1/1 


AD-UCS  C91 
UNCLASSIFIED 


THE  ESTIMATION  OF  QEOPOTENTIRLS  BV  HAV  OF  GEOPHYSICAL 
INVERSE  THEORVCU)  HASSACHUSETTS  INST  OF  TECH  LEXINGTON 
LINCOLN  LAB  H  T  LANE  ET  AL.  27  JAN  86  TR-73S 
ESD-TR-83-279  F19628-83-C-00e2  F/O  8/3  NL 


AD-A165  691 


B0WTMM79 


TedmiaU  Report 
735 


The  Estimation  of  Geopotentials 
by  Way  of  Geophysical  Inverse  Theory 


M.T.  Lmne 
ELM.  Gapo«chkin 


27  January  1986 


Lincoln  Laboratory 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
Lexiwron, 


Prep«rffd  for  ihr  Drpartmrnl  of  the  Air  Forte 
undrr  Eiectronie  Syelrmi  DivUioo  <.ontnici  F19628-a5-C*<K)02, 


Approved  for  public  release;  distribution  unlimited. 


TIm  work  ntponeti  ia  t|iU  tftKumeat  wj|«  >erlora[i*4-iH  Unmin  LatenMr;,  • 
«aater  te»  fMeareh  operated  ky  Matwrkmtu  laaittWa  of  TeahMlagy^  adtk  tke 
(Upport  of  tke  Oepanmeat  of  the  Air  Forte  aador  GaatfaciFlMII  n  C  WK. 

TbU  report  may  he  reproduced  ioeatlafy  aeod*  of  U^S.  Gotauinatal  apeadea. 


The  «iew«  and  conriualon*  conlatned  in  thia  docuateni  are  thote  of  the  contractor’ 
and  ahould  not  be  interpreted  a*  nercaaarily  repeeteatiog  the  official  pollciet. 
either  expreaaed  or  implied,  of  the  CnHed  $late<  CovernmeoL 


The  ESD  Public  Affain  OffichhaK  reviewed  thir  report,  and 
it  i»  releacable  to  the  Nellonai  Technical  Information 
Servi<-e.  where  it  will  be  available  to  the  gea^al  piiMic, 
including  foreign  itationaU. 


Thi.i  technical  report  ha«  been  reviewed  and  i«  approved  for  publicatioa. 
FOR  THECOMMANDEB 

ThomikH  J.  AlpiTi.  Major,  I  SAK 

(Jhief,  ESI)  Linrohi  lr«l>or«lor\  rrojcrl  Offne 


Non  Uncoln  Recipients 

PLEASE  DO  NOT  RETURN 

Permission  is  given  to  destroy  this  document 
vdien  it  is  no  longer  needed 


MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
LINCOLN  LABORATORY 


THE  ESTIMATION  OF  GEOPOTENTIALS 
BY  WAY  OF  GEOPHYSICAL  INVERSE  THEORY 


M.T.LANE 
E.M.  GAPOSCHKIN 

Group  9J 


TECHNICAL  REPORT  73S 

27  JANUARY  1986 


DTK 

SELECTI 
MAR2  4ig8 

B 


Approved  for  public  release;  distribution  unlimited. 


LEXINGTON 


MASSACHUSETTS 


ABSTRACT 

Satellite  to  Satellite  Tracking  data  (SST)  can  be  used  to  measure  the  geopotential 
at  the  satellite  altitude.  This  measurement  can  be  used  to  estimate  the  Earth's  gravity 
field  at  the  Earth's  surface,  the  so-called  “inverse  problem.”  Geophysical  inverse 
theory  is  applied  to  this  inverse  problem,  and  numerical  methods  are  developed 
and  tested.  Geophysical  inverse  theory  is  used  to  map  the  geopotential  from  the 
satellite  altitude  to  a  lower  surface.  Two  configurations  are  explored  and  the  geopotential 
in  a  local  network  is  recovered  with  less  than  4%  error. 
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THE  ESTIMATION  OF  GEOPOTENTIALS 
BY  WAY  OF  GEOPHYSICAL  INVERSE  THEORY 


I.  INTRODUCTION 

In  general,  inverse  theory  begins  with  an  integral  equation 

f(x)  =  J  g(x.y)  u(y)  dy  (I) 

where  the  kernel  g(x,y)  is  given  a  priori,  and  x  and  y  may  be  vectors.  When  u(y)  is  known,  f(x) 
is  determined  by  simple  integration  methods  and  is  known  as  the  “forward  problem.”  Conversely, 
when  f(x)  is  given,  one  seeks  to  determine  u(y)  which  is  known  as  the  “inverse  problem.”  In  our 
case.  f(x)  is  not  known  everywhere,  and  therefore  we  can  only  hope  to  obtain  a  less  complete 
description  of  u(y).  Nevertheless,  one  can  develop  a  mathematical  structure  that  allows  the  max¬ 
imum  information  (in  some  sense)  to  be  obtained  from  a  given  set  of  data  and  a  quantification 
of  the  quality  of  the  result.  Such  a  formalist  i  has  been  developed  by  seismologists  in  solid  earth 
geophysics,  and  is  known  as  “geophysical  inverse  theory.” 

Geophysical  inverse  theory  has  made  a  profound  impact  on  the  study  of  the  Earth’s  interior. 
The  reader  is  encouraged  to  see  Reference  I  or  2  for  an  overview  of  the  theory  and  how  it  has 
been  useful  in  the  geophvsical  context.  A  motivation  for  geophysical  inverse  theory  is  that  much 
of  our  knowledge  of  the  Earth's  interior  is  based  on  observations  made  at  the  surface.  This  is 
opposite  to  the  forward  problem  in  which  a  valid  mathematical  model  is  constructed  from  which 
the  behavior  of  the  interior  can  be  calculated  directly.  Inverse  theory  is  usually  not  capable  of 
extracting  information  that  is  not  intrinsically  contained  in  the  data,  and  therefore,  the  theory  is 
extremely  sensitive  to  corruption  in  the  data  and  the  model  a.ssumptions  from  which  the  data  has 
been  taken. 

The  origin  of  our  inverse  problem  is  as  follows.  For  a  number  of  important  geophysical 
investigations  from  solid  earth  geophysics  to  physical  oceanography,  one  needs  to  know  the  grav¬ 
ity  field  of  the  Earth,  or  geopotential,  with  a  greater  degree  of  accuracy  and  resolution  than  has 
been  achieved  to  date.-^  The  technique  of  Satellite  to  Satellite  Tracking  (SST)  can  provide  mea¬ 
surements  of  the  geopotential  at  the  satellite  height  which  will  be  described  below'.  This  satellite 
data  can  provide  a  global  map,  at  satellite  altitude,  of  the  geopotential.  With  such  a  map. 
methods  are  needed  to  determine  the  potential  at  (or  near)  the  Earth's  surface. 

One  satellite  mission  under  consideration  is  for  two  satellites  in  circular  polar  orbits,  in  the 
same  plane,  at  the  same  altitude  of  160  km,  separated  by  hundreds  of  kilometers,  tracking  each 
other  with  an  accuracy  of  10'^  m  sec.  To  be  workable  it  is  assumed  that  both  satellites  must  be 
equipped  with  surface  force  compensation  (so-called  drag  free)  devices.  An  alternative  approach 
also  under  consideration  would  be  the  use  of  a  gravity  gradiometer  in  a  similar  orbit.  In  either 
case,  the  gravity  field  at  the  Earth's  surface  is  ultimately  needed,  and  the  inverse  methods  de¬ 
veloped  here  would  be  applicable. 


This  report  will  discuss  the  application  of  geophysical  inverse  theory  to  computing  the  geo¬ 
potential  of  points  above  the  Earth  with  respect  to  mass  points  located  on  the  Earth’s  surface. 

The  observed  data  in  this  context  will  be  actual  geopotentials  of  points  located  at  a  greater  dis¬ 
tance  from  the  Earth  than  the  points  at  which  the  geopotentials  are  to  be  estimated.  This  is 
known  as  the  “downward  continuation  problem.”  We  begin  by  describing  the  Earth’s  potential 
(geopotential)  at  any  point  (r,d>,X)  exterior  to  any  mass  as 

V(r,d>,X)  =  GM/r  +  T(r,d.,X)  (2) 

Here  0  ^  ^  jr  is  the  spherical  angle  of  latitude  (assuming  <t>  =  tt/  2  is  the  equator),  0  <  X  ^  27r 

is  the  angle  of  longitude,  G  is  Newton’s  gravitational  constant,  and  M  is  the  mass  of  the  Earth. 
The  central  force  term,  GM/r,  is  chosen  so  that  the  average  of  T(r,d),X),  the  anomalous  potential, 
over  any  external  sphere  is  zero.  Now,  as  general  practice  in  physical  geodesy,  we  concern  our¬ 
selves  with  T. 

From  potential  theory,'*  if  one  knows  the  potential  everywhere  on  a  sphere,  then  the  poten¬ 
tial  at  any  point  outside  that  sphere  can  be  obtained  from  Poisson’s  Integral  Formula  as 

R(r2-R2)  r27r  ptr  T(R,d,'.X')  .  . . 

T(r.d),X)=  — -  J  J  — - - sin«^)'dd)dX  (3) 


where 


p2  =  r2  +  R2  -  2Rr  cos(i/>) 


cos  =  cos  d)  cos  <t>'  +  sin  <t>  sin  d>'  cos  (X  -  X') 

This  is  our  forward  problem. 

From  celestial  mechanics,  one  can  compute  the  perturbation  in  velocity  along  the  direction 
of  motion  (v)  due  to  the  anomalous  potential  (T)  as 

V  =  —  -  2n6r 
na 


where  n  is  the  satellite  mean  motion,  a  is  the  satellite  semi-major  axis,  and  6r  is  a  complicated 
integral  which  depends  on  T  among  other  things.  It  is  known  that  T  na  >  2n5r.  Therefore,  we 
begin  with 

T  =  na(v  +  2n6rj,)  ,  (4) 

where  (Sr^,  is  an  approximation  of  dr  based  on  partial  knowledge  of  T.  In  any  event,  continuous 
observation  of  T  can  provide  a  map  of  T.  and  associated  variances.  We  now  wish  to  use  T(r.d),X) 
to  estimate  T(R.5.X)  on  (or  near)  the  Earth’s  surface. 

We  expect  that  the  density  of  data  and  horizontal  resolution  desired  preclude  inverting  the 
whole  problem  in  one  step.  For  example,  for  one  degree  square  resolution,  more  than  42,000 


parameters  would  be  necessary.  If  orthogonal  functions,  such  as  spherical  harmonics,  were  used, 
more  than  64,000  coefficients  would  have  to  be  determined.  One  motivation  for  using  inverse 
theory  is  to  treat  the  inverse  problem  locally.  That  is,  to  estimate  the  geopotential  at  (or  near) 
the  Earth's  surface,  from  nearby  points  on  an  exterior  sphere. 

It  is  clear  that  our  inverse  problem  of  determining  Tn(y)  =  T(R,d>,\)  from  observed  data 
Tj(y)  =  T(r,<6,X)  is  linear  and  can  be  thus  attacked  by  any  number  of  linear  inverse  methods.  Two 
methods  in  particular  which  have  acquired  familiarity  arc  the  Backus-Gilbert  formulation^.*  and 
the  spectral  expansion  method.^.^.’  The  Backus-Gilbert  approach  requires  an  approximation  of  an 
appropriate  delta  function.  This  approach  is  not  so  practical  in  our  situation  for  the  construction 
of  a  reliable  and  implementable  solution  since,  for  each  point  y  =  (d>,X)  at  which  the  geopotential 
Tr  is  desired,  one  must  approximate  a  delta  function  about  y.  Their  method  does  provide,  how¬ 
ever,  a  means  for  evaluating  the  significance  of  the  solution  (see  Reference  1,  pp.  43-44).  The 
spectral  expansion  method  is  much  more  suited  to  the  construction  of  a  solution  which  can  be 
subsequently  evaluated  at  several  points,  y,  and  we  will  pursue  this  method  in  our  inverse 
problem. 

Most  solutions  by  inverse  methods  are  not  unique,  and  it  is  necessary  to  briefly  discuss  this 
problem  in  our  context.  The  question  of  uniqueness  in  our  formulation  is  equivalent  to  the  ques¬ 
tion  about  the  nature  of  the  “annihilalor  of  the  kernel.”  Precisely,  the  annihilator.  A,  is  defined 
by 

A  =  {b(y)  I  J*  g(x,y)  b(y)  dy  =  O} 

Sr 

where  Sr  is  the  sphere  of  radius  R.  If  A  is  empty,  then  the  solution  T^fy)  is  unique.  We  cannot 
count  on  this  being  the  case,  however,  and  in  fact,  A  turns  out  to  be  infinite  dimensional.  We 
must  proceed  without  necessarily  having  uniqueness,  however,  and  simply  be  aware  of  this  prob¬ 
lem  which  could  offset  the  data.  One  possible  consequence  of  nonuniqueness  is  that  even  the  sign 
of  the  estimated  potential  may  be  wrong. 

This  report  contains  five  sections,  including  the  introduction.  In  Section  11  we  will  outline 
the  spectral  expansion  method  and  discuss  some  philosophies  concerning  its  use.  Section  111  con¬ 
tains  information  regarding  the  physical  structure  of  the  experiments  performed  using  the  spectral 
expansion  technique.  Section  IV  follows  with  a  detailed  discussion  on  the  computational  aspects 
of  the  experiments,  and  we  conclude  with  Section  V,  which  outlines  the  numerical  results  and 
some  conclusions  that  can  be  made. 


V 


II.  THE  SPECTRAL  EXPANSION  METHOD 


In  this  section,  we  will  discuss  the  spectral  expansion  method  and  introduce  the  tools  that 
will  be  needed  for  construction  of  the  solution  Tj^.  We  will  follow  Parker's  outline  (Reference  I, 
pp.  46-49).  Suppose  we  have  a  finite  number,  N.  of  observed  data,  each  having  the  form  in  (3). 
Then  we  can  write 

Tr(x,)  =  J  g(Xi,y)  TrIv)  dy  i  =  I . N 

Sr 

or 


gi(y)  TR(y)  dy 


i  =  I . N 


(5) 


If  we  need  to  treat  corrupt  data,  then  we  can  assume  a  distribution  for  the  error  and  include 
the  "noisy"  data  from  the  beginning.  We  consider  (5)  weighted  by  the  inverse  of  the  standard 
error  of  each  measurement,  a^: 


gi(y) 


I  R(y)  dy 


which  we  will  write  as 

=  X  Tfi(y)  dy 

Sr 


i  =  I . N 

i  =  I . N 


(6) 


Now  T'j  is  dimensionless  with  unit  variance. 

We  proceed  by  introducing  a  matrix  P  with  elements  Pjj  formed  by  integrating  the  double 
kernel 


•'ij  =  X  Slly)  gj(y)  dy  (7) 

Sr 

r  is  non-negative  definite  and  symmetric,  and  thus,  it  can  be  diagonalized  by  an  orthogonal 
matri.x  0, 

O'rO  =  .\  =  diag(A, . (8 

with  A|  >  Ai  ^  ^  ^  0.  Zero  eigenvalues  may  occur  if  some  gj  is  a  multiple  of  some  gj 

with  i  ^  j,  but  we  will  disregard  this  possibility  temporarily.  The  set  of  eigenvalues  is  called  the 
"spectrum”  of  the  problem  defined  by  (6). 

Consider  the  functions 

N 

'/'i(y)  ^  '  1  d|,g;(y) 


(9 


It  is  easy  to  verify  that  the  t/>|’s  are  orthonormai  functions;  that  is 


J  '/'j(y)  dy  =  Sjj 
Sr 

where  6jj  is  the  Kronecker  delta  function.  Thus  we  may  expand  Tr  in  terms  of  these  functions 
N 

TR(y)  =  S  ajiAify)  +  iA*(y)  (10 

i=l 

where  i/r*  belongs  to  A,  the  annihilator  of  the  kernel,  and 

j*  ^*(y)  ^>(y)  dy  =  0 
Sr 

for  each  i.  The  coefficients,  a;,  of  this  expansion  must  be 

N 

ai  =  S  '('i(y)  TR(y)  dy  =  X-'  ^  X  0|j  /  gjfy)  T^fy)  dy 
Sr  j=|  Sr 

N 

=  A-'2^0ji-rj  (11 

j=l 

Parker  (Reference  1,  p.  47)  notes  that  the  standard  error  of  each  coefficient,  aj,  is  X‘*  -  and 
that  each  aj  is  statistically  independent.  Thus  the  coefficients  of  the  orthonormal  functions  lAj 
increase  in  uncertainty,  and  the  uncertainty  for  «/>•(>’)  is  total.  If  any  zero  eigenvalues  appear,  the 
functions  associated  with  those  eigenvalues  should  be  included  in  the  function  0*  undetermined 
by  the  data  because  those  functions  have  coefficients  which  are  complete  in  uncertainty.  The 
functions  also  become  more  oscillatory  as  i  increases,  and  so  the  smoothest  parts  of  Tr  are 
most  accurately  determined  because  the  coefficients  of  the  smoother  functions,  i/»j.  are  more  pre¬ 
cise.  Thus,  when  many  of  the  eigenvalues  are  a  small  fraction  of  X|.  then  we  can  assume  that  we 
have  a  numerically  (and  physically)  poorly  posed  problem,  since  in  this  case  many  of  the  coeffi¬ 
cients  a;  exhibit  large  variances. 

It  is  conceivable  to  try  and  provide  a  smoothed  version  of  the  true  Tr.  and  we  can  do  this 
by  filtering  out  the  functions  whose  coefficients  have  large  uncertainties.  Thus,  we  have  a  trun¬ 
cated  version  of  ( 10) 

L 

TR(y)  =  ^  a>,(y)  +  !/>*(>)  wiih  I.  N  .  (12) 


Unfortunately,  it  is  not  easy  to  determine  the  proper  choice  of  L  which  maximizes  the  accuracy, 
since  the  truncated  sum  (12)  no  longer  fits  the  data  exactly.  Parker  (Reference  1,  pp.  48-49) 
introduces  the  squared  two-norm  misfit  of  the  data  as  a  possible  tool  to  study  the  effects  of  trun¬ 
cation.  Define 


where 


^  S  (f;  -  t;)2 

i=l 


T,'  =  S  gi(y)  TR(y)  dy  for  each  i 
Sd 


Using  (7)  and  (8),  we  can  derive  a  more  practical  expression  for  x^-  From  (13),  we  have 


i=l' 

U=L+1 

/  N 

--  x{ 

s 

i=l' 

^ j=L+l 

N 

/  N 

= 

X 

i=l' 

j=L+l 

N 

!  N 

= 

S 

i=l ' 

V  j=L+l 

S  ajJ  gi(y) '/'j(y)  dy ) 

L+1  Sr  / 

N  N  ' 

X  SajX:'  \j  r  gi(y)gk(y)dy 


N  N  N 

=  S  X  SajaKAJ  2x1/2oTo,j^ 

j=L+l  k=L+l  i=l 

N  N 

=  S  S  aja^x;  2xjl  2(0T0)., 

j=L+l  k=L+l 

N 

=  X  a/Xj  (14) 

j=LH 

Parker  (Reference  1,  p.  49)  suggests  choosing  L  s  that  equals  its  expected  value,  or  so  that 
P(x^)  -  b.5.  Another  suggestion  is  to  choose  L  so  that  x^  is  close  to  the  number  of  independent 
data.  One  must  proceed  with  caution,  however,  since  the  proper  choice  of  L  is  intimately  related 
to  the  model  of  the  underlying  system,  as  the  reader  will  discover  in  Section  V  (where  we  outline 
the  results  of  our  experiments). 
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III.  THE  APPLICATION  OF  THE  SPECTRAL  EXPANSION  METHOD 


Two  experiments  were  conducted.  For  each  experiment,  a  model  problem  was  developed;  the 
first  to  accommodate  a  global  setting,  and  the  second  to  accommodate  a  local  setting. 


In  the  first  configuration  (representing  the  global  setting),  twenty  mass  points  were  spaced 
evenly  on  the  Earth’s  equator.  Each  point  was  assigned  a  unit  mass  with  alternating  sign.  Thus, 
the  mean  mass  anomaly  is  zero.  Then,  N  sample  points  were  chosen  randomly  at  a  constant 
radius  above  the  Earth  within  a  10%  latitude  band  of  the  equator  =  njl).  The  geopotential 
of  each  sample  point  was  then  computed  with  respect  to  the  twenty  mass  points.  Thus,  the 
observed  data  for  the  i^^  sample  point  is 


Ti  =  G 


20 


2 

j=l 


(15) 


where  Rjj  is  the  distance  between  the  (i^l^)  sample  point  and  the  mass  point.  G  is  Newton’s 
gravitational  constant.  One  unit  of  distance  was  chosen  to  be  the  radius  of  the  Earth,  and  there 
G  =  1.5362E-06  units^  sec'^.  A  rough  sketch  of  configuration  1  is  shown  in  Figure  1. 


Figure  I.  Global  picture  of  configuration  I  with  mass  point  structure. 

In  the  second  configuration,  representing  the  local  setting,  we  chose  a  more  condensed  struc¬ 
ture.  One  mass  point  was  placed  at  the  intersection  of  the  meridian  and  the  equator,  and  a  grid 
of  24  other  mass  points  was  placed  around  this  “origin”  with  the  points  spaced  at  1°  intervals. 
Each  point  was  assumed  to  have  unit  mass  with  alternating  sign  as  in  the  first  configuration.  The 
structure  is  shown  in  Figure  2  with  the  origin  circled. 


.■TWjvr- 


EQUATOR  ^  =  7r/2 


Now  121  sample  points  were  set  up  in  a  similar  grid  (11  in  each  row  and  column)  at  a  con¬ 
stant  radius  above  the  Earth  directly  above  the  grid  of  the  2S  mass  points.  Hence  the  origin  of 
the  sample  point  grid  is  also  at  0  longitude  and  vjl  latitude.  The  geopotential  of  each  sample 
point  was  then  computed  with  respect  to  the  25  mass  points  similar  to  (15). 

Since  the  mean  mass  anomaly  of  the  second  configuration  is  not  zero,  we  created  an  alter¬ 
nate  configuration.  We  removed  the  mass  at  the  origin  in  the  grid  of  the  25  mass  points  and  per¬ 
formed  an  alternate  experiment.  In  order  for  the  grid  to  be  “isotropic”  (invariant  under  rotation 
and  orientation),  the  mass  points  were  all  assigned  the  exact  sign  as  in  the  configuration  using  all 
25  of  the  mass  points  (see  Figure  3).  Sample  points  were  set  up  in  the  same  way  as  before. 

With  each  configuration,  at  a  radius,  R,  strictly  between  the  Earth’s  radius  and  the  smallest 
radius  at  which  a  sample  point  was  taken,  the  spectral  expansion  method  is  called  to  yield  the 
geopotential  Tg.  We  ignored  in  (10)  and  (12),  since  it  is  undetermined  by  the  data. 


IV.  NOTES  ON  COMPUTATIONS 

Before  we  analyze  the  results  obtained  for  the  experiments  performed,  it  is  necessary  to 
understand  exactly  how  the  computations  were  derived.  Indeed,  since  the  computations  were 
complicated  and  involved  large  matrices,  one  may  look  upon  the  error  in  the  results  with 
some  caution  rather  than  if  the  computations  were  straightforward.  Once  the  matrix  P  is  found, 
the  spectral  expansion  method  required  the  computation  of  eigenvalues,  eigenvectors,  and  several 
finite  sums. 

The  routine  used  to  compute  the  eigenvalues  and  eigenvectors  was  the  well-known  QL 
algorithm,  a  close  relative  of  the  QR  algorithm.  Here  L  is  tridiagonal,  which  was  converted  from 
the  original  matrix  P  by  Householder’s  method  (see  Reference  8).  This  method  was  chosen  for  its 
speed  and  economy  of  use  of  storage  space.  In  the  context  of  large  matrices  where  all  the 
eigenvalues  and  eigenvectors  are  needed,  we  found  this  method  to  be  more  accurate,  much  faster, 
and  requiring  less  storage  space  than  other  methods  such  as  bisection  or  the  Jacobi  method. 

It  remains  to  discuss  the  computation  of  P,  (7),  and  we  will  give  a  detailed  account.  We 
write  (7)  explicitly  as 

rp2ir 

=  J  gfrj.  <l>v  *■<  r.  d>,  X)  R2  sin  </>  d\dd> 


-■f-i  - - 

0  0  I47r  RpfJL47r  Rpj, 


dXdi^ 


p?=  r?~  R2  -  2ri  R  cosilrjy  and  p?=  r?+  R2  -  2rj  R  cosi^jy 

where  ^jy  is  the  spherical  angle  between  (i^j,  Xj)  and  y  =  (</>,  X).  We  may  rewrite  cos  i/rjy  using  the 
following  well-known  law  of  spherical  trigonometry; 

cos  i/fjy  =  cos  i/fjj  cos  i/fiy  +  sin  i/»ij  sin  «/rjy  cos  a  ,  (17) 

where  a  is  the  angle  between  the  arcs  and  flf^y  (see  Figure  4).  Now  as  0  and  X  range  over  the 
whole  sphere,  a  and  i/fjy  will  also  range  over  the  whole  sphere  with  0  ^  a  ^  2n-  and  0  ^  i/fjy  ^  tt. 
Hence,  integrating  the  double  kernel  with  respect  to  a  and  i/f^y  will  yield  an  equivalent  integral, 
and  we  can  thus  write  (16)  as 


rr 

■rf-R2- 

tj- R2 

sin  1 

0  0 

.47r  Rpf. 

.4n-  Rpj^. 

=  J  J  - :  - -sini/Tj  dad</»i 

0  0  l47rRpfJl4a-Rp^ 

We  will  proceed  to  simplify  this  expression,  but  first  we  rename  constants  as  follows: 
(rf-R2)(r?-R2) 

c, - c,  =  ,f,R^c,=  r,*R^.c,  =  ^2r,R, 

Cj  =  -2rjR  ,  c^  =cos  tj/-  and  Cy  =  sin 


l«Pj.  Aj) 


v  =  «^.  X) 


Figure  4.  Spherical  triangular  relationship  among  points. 

We  will  also  write  tit  -  tli^y.  Thus,  we  write  (18)  as 
>  C|  sin  til  dadtii 

:os  til)  (C3  +  C5  (c^  cos  tit  +  c-j  sin  ili  cos  a))]  3/2 

COS0,  a2(iA)  =  C3  +  C5C6  cos  lit,  and  a){tli)  =  c^ci  siniA,  we  have 

C|  sint^  dadtit 

_  _ 


^(C2a2  (tit)  +a((i/f)  a2(i/r))  (1  +  C0Stt)j3/2 


To  simplify  (17)  further,  we  set  b|(i/>)  =  C|  sin  i/r/ [C2a2 (t(r)  +  a|(i/r)  a2(i/»)]3/2  and 
b2(tli)  =  a3(«/f)/a2(i/f).  Then, 


rij=  r  b|(*)  r 


0  0  (1  +  b2(i/r)  cos  a)3  2 

We  leave  (21)  temporarily  and  examine  the  nature  of  b2(iA)-  Indeed,  we  will  show 
that-l<b2(tl»)^0  for  all  OKtli^jr.  Recall 


h2(tli)  = 


-2rj  R  sin  ili^^  sin  tit 


r?+  R2  -  2rjR  cos  t|l^j  cos  ili 

Since  0  ^  tit  ^  n,  the  numerator  of  (22)  is  always  negative  (since  rj,  r  >  0).  Moreover,  -1  ^ 
cosi(fij  cosi/r^l  and  rj>r  imply  that  the  denominator  of  (22)  is  always  positive  for 
since  r?+  R2  >  2rjR.  Thus,  bj  (tit)  is  always  negative  for  any  0^  tit  Now,  since  -1  ^ 
cos(tli,,-ilt)^\,  it  follows  that 


67501 -N 


2rjR  cos  -  lA)  =  2rjR[cos  cos  >1/  +  sin  i^ry  sin  ^^^]<  i|  +  R^  , 
which  implies  that 

i^+  R2  -  2rjR  cos  cos  tl/  >  2rjR  sin  sin  ^  for  all  ^^y,  ^ 

Hence  b2i^)>  -1  as  claimed. 

We  ignore  ^  for  the  moment,  and  we  set  a  =  -b2itlt)  (so  that  0<  a<  1). 

We  wish  to  simplify  further: 

r2n  da  _  _ do _ 

0  (1- a  coso)3/2  Q  [1  -  a(l  -  2  sin2^)]3/2 

f"”  2du  r  a, 

=  J  - - — —  [where  u  = — 1 

*0  ((1  -  a)  +  2a  sin2u)3/2  2 


[where  k  =  •^] 


4 

pjr/2 

du 

'  (1  -a)3/2 

(1  +  2a  sin2  u/(l  -a))3/2 

4 

pn-/2 

du 

"  (1  -a)3/2 

(1  +  k(l  -  cos 2  u])2/2 

-4 

f  ® 

dy 

(1  -  a)2/2 

■  ir/2 

((1  +  k  -  k  sin2y)3/2 

_ 4 _  _ 

(1  -  a)^/2  (1  +  |c)3/2  (1 -Tc^  sin^y)^/^ 


[where  y  =  y  u] 


[where  ^ 

k  +  1  a  +  r 


Now  (24)  is  an  elliptic  integral  of  the  third  kind,  and  if  we  use  formula  17.7.24  in  Reference  9, 
we  are  able  to  reduce  it  to  an  elliptic  integral  of  the  second  kind  which  can  be  approximated 
numerically  to  any  degree  of  accuracy  by  using  well-known  techniques.  Indeed.  (24)  can  be  writ¬ 
ten  as 


4  pn  _  _ 

— Tjn  ^  J  V  I  -  k^  sin^y  dy  [where  sin  ^7  =  k^] 

Since  sec^y  =  (I  ♦  a)/(l  -  a),  we  take  as  our  final  form  of  (23) 


V(«  = - i™/’ 

( I  -  a)  V  I  +  a  0 


2  /  2a  .  ,  ^ 

/I  -  ; -  sin-y  dy 

v  I  +  a 


4  ^  /  2b2(dr) 

- ,,  .  .  -  J  /I  +  -  sin^y  dy 

(I  +  bjiti/)  \/ 1  -  b2  (ill)  0  V  I  b2  (i/») 


where  I  <  2b2(i/f)/(l  b2(i/»))  ^  0. 


There  is  a  vast  amount  of  literature  on  elliptic  integrals,  and  computation  to  any  degree  of 
accuracy  desired  is  possible  by  a  variety  of  techniques.  We  will  outline  below  two  methods  for 
computation  of  Y(i^),  for  a  given  t/»,  which  are  particularly  good;  a  (pseudo)-hypergeometric  series 
expansion  and  an  expansion  in  terms  of  theta  functions.  We  will  concentrate  on  the  following 
general  form  of 

E(k)  =  x'^yr  -  sin^  y  dy  (27) 

0 

where,  in  our  instance,  =  -2b2(iA)/  (l  -  b2(tl»))  for  a  given  Notice  that  0  ^  k^  <  1. 

We  begin  with  Grobner  and  Hofreiter  (Reference  10,  formula  2b,  2c,  p.59) 


E(k)=  '[^(-k2)  f2„(^) 


n=0  n 
where  f^lTr  2)  =  7r;2  and 
2 

2)-  J 

^  0 

Thus, 


hn  =  J  sm2ny  dy  =  f2„.2  (:;) 


2n 


f2n(l) 


(I;  2;  n)" 


2'  2"*'  n! 

where  (1;  2;  n)  =  1(1  +  2)  (1  +  4)  (1  +  6)  .  .  .  (1  +  2(n  1)).  Now 


('  ^)  (1  2)"  =  [(1)  (  -  2)  (1  -  4)  ...  (1  2(n  -  1))] 

n .  n' 


and  so  if  we  substitute  (29)  and  (30)  into  (28).  we  have 


n=0  22'’*'  (n!)2 


[2(n  l)]2)] 


(28) 


(29; 


(30) 


(3i: 


Equation  (31)  is  the  (pseudo)-hypergeometric  series  representation  of  (27).  The  ratio  test  of  (31) 
yields  k2,  and  since  0  ^  k2  <  1,  it  does  follow  that  (31)  will  always  converge.  If  k2  is  quite  small, 
then  (31)  will  be  sufficient  to  give  any  degree  of  accuracy  desired  with  only  a  few  terms  in  the 
expansion,  but  as  k2  tends  more  toward  its  upper  bound,  the  scries  (31)  converges  slower,  mak¬ 
ing  it  inappropriate  for  all  possible  situations.  This  leads  us  to  a  more  universal  representation  of 
(27);  by  using  theta  functions. 


We  begin  by  listing  the  four  most  common  theta  functions  (Reference  1 1,  p.  464).  For  /  a 
complex  number  and  |q|  <  I,  define 


Jl|(z,q)  =  2q''^  sinz  -  sin3z  +  2q"/'»  sinSz  .  .  . 


=  2  S  (-I)"  q<"*'  2)2  sin[(2n  +  I)z] 
n=0 

Jl2(z.q)  =  2q'  ^  cosz  +  2q^  ^  cos3z  +  2q25  4  cosSz  .  .  . 


=  2  S  q<"*'  2)^  cos  [{2n  +  1)  z] 
n=0 

il3(z,q)  =  1  +  2q  cos2z  +  2q^  cos4z  +  2q^  cos6z  .  .  . 


=  1  +  2  S  q"^  cos(2nz) 
n=l 

and 

8-4(z.q)  =  1  -  2q  cos2z  +  2q'‘  cos4z  -  2q^  cos6z  +  .  .  . 


=  I  +  2  S  (  1)"  q""  cos(2nz) 
n=l 

Denote  *"<1  <i4(0,q)  by  £2.  *-3'  and  ?,4.  respectively,  and  write  H'l  and  l\'  to  denote 

the  first  and  third  derivatives,  respectively,  of  il|(z,q)  with  respect  to  z.  evaluated  at  z  =  0.  Whit¬ 
taker  and  Watson  (Reference  11,  p.  467)  note 

and  so  if  we  set  k  =  complement  of  k.  k  (k2  +  k2  =  |),  satisfies,  v/T  =  ll4/!2.3.  We 

form  2e  =  (1  -  \/T)/(l  +  x/lT)  and  note  that  0^  t  <  1/2  whenever  0^  k  <  1.  We  follow  Refer¬ 
ence  11,  p.  486  by  using  t  to  represent  q: 


q  =  e  +  2t5  +  I5t9  +  i50el3  +  0(t'’)  (3: 

which  converges  very  quickly  (usually  after  only  the  first  three  or  four  terms).  Now,  if  we  define 
K  =  1, 2n-Jl^,  then  it  can  be  shown  (Reference  1 1,  pp.  499-500),  that 


which  is  an  elliptic  integral  of  the  first  kind.  Finally,  it  can  be  shown  that  (27)  can  be  written  as 


Now  we  are  given  k^,  and  so  we  can  use  (33)  to  find  an  accurate  value  for  q,  and  then  proceed 
to  find  I'l',  £3,  and  i\  which  can  be  used  to  solve  for  E(k)  by  using  (35).  This  is  the  best  formula¬ 
tion  we  have  found  to  yield  (27)  as  accurately  as  possible  given  all  values  of  0  ^  ^  tt. 

We  continue  with  our  development  of  computing  Tjj  in  its  final  form; 

b|(i/f)  Y(i/r)  di/r  (36) 

0 

We  evaluated  (36)  numerically  using  cautious  Romberg  extrapolation  to  find  an  acceptable 
estimate  of  the  integral  on  various  subintervals  suitably  chosen  over  the  given  interval  0^  i/r  ^  rr. 
We  simultaneously  solved  for  Y(iA)  by  way  of  either  (31)  or  (35),  depending  upon  t/r.  For  both  of 
these  integrations  we  used  subroutines  developed  by  the  Internationally  Mathematical  and  Statis¬ 
tical  Library  (IMSL). 

Instead  of  using  the  Romberg  method  on  the  outer  integral,  we  could  have  used  the  Gaus¬ 
sian  quadrature  technique  on  the  interval  [0,  tt].  Ihis  method  was  tested,  and  it  was  found  that 
this  technique  took  much  longer  than  the  former  method  .since  the  integrand  is  not  always  well 
behaved.  The  Romberg  method  is  often  capable  of  handling  jump  discontinuities.  Whatever 
method  used,  one  should  keep  in  mind  that  approximating  (36)  is  more  complicated  and  takes 
longer  as  either  r^  or  rj  approaches  R,  the  radius  of  the  sphere  of  integration.  This  may  effect  the 
results  of  the  spectral  expansion  method  when  sample  points  are  chosen  close  to  the  points  where 
3  is  to  be  estimated. 

One  final  note  on  computation  of  the  matrix  F.  If  we  fix  r,  =  rj  =  c,  a  constant  lor  all  i  and 
j,  then  it  is  clear  that  each  element  I'n  will  only  depend  upon  the  angle  between  the  points 
P,(c.  d)|,  A|)  and  Pi(c.  <t>y  A,).  There  docs  not  seem  to  be  a  closed  form  solution  for  (36)  in  this 
situation,  but  one  method  of  computing  L  more  rapidly  is  to  calculate  I’lj  for  a  fraction  of  the 
total  number  of  matrix  elements,  lit  a  curse  through  those  values  corresponding  to  the  angle 
between  P,  and  P^  and  interpolate  to  tigure  the  elements  I’jj  that  were  not  directly  computed,  \^e 
did  not  do  this  in  our  experiments,  but  this  could  reduce  the  computation  time  considerably  it  a 
large  number  of  sample  points  was  insolved. 
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V.  NUMERICAL  RESULTS 

In  each  experiment,  we  had  the  advantage  of  knowing  the  true  value  of  the  predicted  Tj^fy) 
by  simply  computing  the  geopotential  with  respect  to  the  mass  points  in  the  same  way  the 
observed  data  were  gathered  as  in  (IS).  This  is,  of  course,  not  always  possible  in  the  situations 
discussed  in  the  introduction,  but  the  models  in  our  experiments  were  specifically  chosen  so  that 
we  could  test  and  understand  the  relevance  of  inverse  theory  in  the  context  of  computing  geopo¬ 
tentials.  Hence,  we  are  able  to  compute  a  measure  of  the  performance  of  the  inverse  theory 
which  is  entirely  analogous  to  the  Backus-Gilbert  method  (Reference  1,  p.  44).  For  each  test 
point  (R,  y^),  after  computing  the  predicted  potential  T^^fyi)  and  the  actual  potential  Tg^fyj),  we 
computed  the  standard  deviation  divided  by  the  average  value  of  the  actual  potentials  over  the 
range  of  sample  points: 

N 


s  = 


S  (TR(yi)  -  TR(yi))2 
i=l _ 

1  T2(yi) 

i=l 


(37) 


Considering  integral,  eigenvalue,  and  eigenvector  approximations,  we  will  be  satisfied  if  our 
results  yield  0  <  S  ^  0.2,  and  we  can  conclude  that  not  enough  information  is  contained  in  the 
data  if  S  >  0.2. 

It  should  be  made  clear  that  the  inverse  routine  is  computationally  expensive,  and  so  before 
we  analyze  the  results  from  the  experiments,  we  will  briefly  discuss  the  routines  and  the  time 
involved  in  executing  the  routines.  The  word  size  of  the  program  with  N  =  500  is  close  to  four 
million,  and  so  the  physical  limitations  of  the  computer  limit  the  amount  of  sample  points  one 
can  choose.  Perhaps  a  computer  that  allows  more  storage  space  and  is  faster  could  make  this 
routine  more  viable. 

Our  routine  was  run  on  the  Harris  H800,  and  Table  1  lists  the  time  involved  for  the  execu¬ 
tion  of  various  aspects  of  the  routine  for  configuration  I  with  different  values  of  N,  the  number 
of  sample  points,  and  with  the  outer  sphere  containing  the  sample  points  as  l.l  and  the  sphere  of 
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TABLE 


TABLE  IV 

N 

Largest  Eigenvalue 

Smallest  Eigenvalue 

20 

24.78332 

7.156385 

50 

45.99 

3.00429E-02 

100 

66.99387 

2.31888E-02 

200 

121.4457 

9.078775E-03 

250 

145.5075 

9  0663478E-03 

300 

147.6130 

8.234201 7E-03 

500 

226.612 

1  2824643E-03 

A 


integration  as  1.05.  From  Table  I,  we  can  see  that  the  computation  of  the  eigenvalues  and  eigen¬ 
vectors  is  slightly  less  than  an  n^  problem.  This  is  not  surprising  since  the  computation  of  the 
eigenvalues  and  eigenvectors  in  this  case  is  equivalent  to  the  computation  of  the  inverse  of  f. 

The  computation  of  T  itself  appears  to  be  slightly  less  than  an  n^  problem,  but  the  computation 
of  r  takes  much  longer  than  the  computation  of  the  eigenvalues  and  eigenvectors  of  T.  However, 
the  computation  time  for  T  will  change  if  one  changes  the  sphere  of  integration  and  the  sphere  at 
which  the  sample  points  were  taken  (which  was  held  constant),  or  if  one  changes  the  configura¬ 
tion  of  the  sample  points.  For  example,  the  CPU  time  involved  for  computing  T  in  configuration 
2  with  121  sample  points  with  the  outer  sphere  1.2  and  the  inner  sphere  of  integration  1.1  was  6 
min  13.523  sec,  which  is  less  than  in  configuration  1  with  only  100  sample  points  taken  with  the 
outer  sphere  and  inner  sphere  1.1  and  1.05,  respectively.  The  computations  in  the  rest  of  the  rou¬ 
tine  (coefficients,  orthonormal  functions,  and  tests)  are  linear,  and  so  if  one  needs  a  satisfactory 
upper  bound  for  the  total  CPU  time  involved  in  the  routine  for  the  first  configuration  involving 
similar  parameters  with  N  different,  then  one  should  use  the  total  CPU  time  of  8  min  49.41  sec 
for  N  =  100,  compute  n  =  N/lOO,  and  form  l/3[n  +  n^  +  n^]  529.41  sec. 

In  the  global  setting,  configuration  1,  there  is  a  direct  relationship  between  the  number  of 
sample  points  used  and  the  tendency  of  S  toward  zero;  however,  the  convergence  is  slow.  Con¬ 
sider  Table  II,  which  lists  S  corresponding  to  varying  N. 

With  each  of  the  cases  listed  in  Table  II,  we  used  the  radius  of  the  outer  sphere  and  the 
radius  of  the  inner  sphere  (of  integration)  as  l.l  (637.8135  km  above  the  Earth)  and  1.05 
(318.90675  km  above  the  Earth),  respectively.  This  assumes  that  the  radius  of  the  Earth  is  1  unit 
(6378.135  km). 

S  was  satisfactory  only  when  N  was  500  (the  upper  limit  of  our  program  size).  So  we  must 
consider  the  effects  of  truncation  as  outlined  in  Section  11.  The  advantage  of  knowing  the  actual 
value,  T[^,  enabled  us  to  study  the  effects  of  using  the  truncated  sum  Tg  in  (12)  and  the  rele¬ 
vance  of  (14).  We  were  able  to  find  the  actual  L  which  gave  the  best  value  for  S.  Table  III 
shows  some  values  of  N  together  with  the  corresponding  L  which  yielded  the  smallest  value  of  S. 
The  same  parameters  of  radii  were  used  as  in  Table  II.  It  is  easy  to  see  from  Table  111  that  an 
elementary  way  to  choose  L  in  all  cases  is  not  likely. 

Before  any  conclusions  can  be  drawn  from  Table  111,  it  may  be  helpful  to  consider  the  spec¬ 
trum  of  the  eigenvalues  for  the  various  values  of  N.  Table  IV  lists  for  the  situations  in  Table  II 
the  corresponding  maximum  and  minimum  eigenvalue.  Figures  5  to  8  show  the  spectra  for  the 
situations  listed  in  Table  III,  where  a  continuous  curve  was  drawn  through  the  eigenvalues. 

Notice  that  not  only  does  the  spectrum  widen  as  N  increases,  but  also,  as  one  can  see  in 
Figures  5  to  8.  that  as  N  increases,  the  graph  of  the  spectrum  resembles  more  of  a  scaled  y  =  I 
curve.  This  widening  of  the  spectrum  indicates  that  more  information  is  contained  in  the  data. 
This  suggests  that,  as  N  increases,  the  more  oscillatory  functions  in  the  expansion  (10)  have  a  less 
damaging  effect  on  the  estimation  of  Tg  in  the  global  setting. 


Figure  8.  Eigenvalue  spectrum  for  configuration  I  with  N  =  200. 


This  is  supported  in  Table  ill.  When  N  is  20,  50,  or  100  (where  the  spectrum  is  relatively 
small),  L  was  less  than  ^09c  of  N  and  the  corresponding  values  of  S  were  greatly  improved.  This 
says  that  in  the  global  setting  when  few  points  are  chosen,  the  oscillatory  functions  with  coeffi¬ 
cients  exhibiting  large  uncertainties  have  a  dramatic  effect  upon  the  estimation  of  Tr-  When  N  = 
200,  L  was  more  than  70%  of  N  and  S  was  only  improved  by  a  small  amount.  Thus,  as  N 
increases,  our  inverse  is  becoming  more  stable  in  the  sense  that  the  solution  depends  continuously 
(with  the  two-norm)  on  the  data.  Thus  the  data  do  not  “smooth"  T^  as  much  as  N  increases, 
indicating  a  closer  fit.  This  observation  is  supported  in  Parker  (Reference  I,  p.  42). 

Finally,  we  consider  the  value  of  x~-  Parker  (Reference  1,  p.  49)  conjectures  that  x-  should 
be  chosen  close  to  the  number  of  independent  data.  This  was  not  found  to  help  in  our  experi¬ 
ments:  for  N  =  20,  50,  or  100,  x^  was  much  larger  than  N,  but  when  N  =  200,  x^  dropped  to  less 
than  1%  of  N.  More  study  is  necessary  before  any  information  from  x’  can  be  extracted  to  help 
us  select  the  L  yielding  the  best  S. 

In  configuration  2,  we  held  the  number  of  sample  points  constant  at  121  and  varied  the  radii 
of  the  sample  point  sphere  and  the  sphere  of  integration.  The  performance  of  the  routine  was 
much  better  over  all  than  in  the  global  setting.  Consider  Table  V  listing  two  tests  of  the  experi¬ 
ment;  the  original  configuration  is  compared  with  the  alternate  configuration  (where  the  mean 
mass  anomaly  was  0). 


TABLE  V 

Radii 

km  Above 
Earth 

S  with  25 

Mass  Points  Used 

S  with  24 

Mass  Points  Used 

(1)  r  =  1.1 

R  =  1 .05 

637.8135 

318.90675 

0.15 

0.019 

(2)  r  =  1 .0333 

R=  1.01667 

212.60429 

106.30182 

0.197 

0.04 

Notice  the  sensitivity  of  the  theory  to  the  model;  the  performance  of  the  alternate  configura¬ 
tion  is  certainly  as  accurate  as  can  be  expected  while,  when  25  mass  points  were  used.  S 
approached  its  upper  acceptable  limit.  Ihe  local  configuration  consistently  out  performed  the 
global  configuration,  c\en  when  the  mean  mass  anomaly  was  not  zero.  We  suspect  the  poor  per¬ 
formance  of  the  nonzero  mean  configuration  is  related  to  the  question  of  uniqueness  and  the 
annihilator.  See  the  later  discussion  of  the  Stokes  integral.  I  his  makes  a  strong  argument  for  the 
local  configuration  (especially  the  alternate  configuration)  as  a  setting  for  a  solution  to  the 
downward  continuation  problem.  Not  only  are  the  results  extremely  satisfactory,  but  also  the  size 
of  the  required  program  is  much  smaller  than  the  size  of  the  program  needed  for  a  similar  accu¬ 
racy  in  the  global  setting.  Furthermore,  truncation  of  the  expansion  for  I  ^  in  the  local  setting 
does  little  to  change  the  value  of  S. 


I  his  can  be  explained  by  once  again  considering  the  spectrum  of  the  eigenvalues.  Because  of 
the  condensed  structure  of  the  sample  points,  the  spectrum  is  much  wider  than  in  a  similar  situa¬ 
tion  in  the  global  setting.  For  example,  when  r=  1.1  and  R  =  1.05,  the  eigenvalues  ranged 
between  7X6.6018  and  I.7I98E-04,  and  the  spectrum  is  .shown  in  Figure  9;  when  r=  1.0333  and 
R  =  1.01667,  the  eigenvalues  ranged  betwen  1843.756  and  3.5439,  and  the  spectrum  is  shown  in 
Figure  10.  One  can  see  from  Figures  9  and  10  that  the  curves  outlining  the  spectra  resemble 
more  of  a  scaled  y  =  I  x  curve,  even  though  only  121  sample  points  were  chosen.  Thus  stability 
of  the  local  inversion  is  better  than  in  the  global  setting  with  similar  parameters.  We  can  con¬ 
clude  that  by  using  all  of  the  eigenvalues  in  the  expansion  (10).  we  will  be  obtaining  a  solution 
close  to  the  best  possible  solution  based  on  the  data. 


Figure  9.  Eigenvalue  spectrum  for  configuration  2  (N  •  121)  with  r  =  /./  and  R  -  1.05. 


TABLE  VI 


N 

Mean  of  Tp 

Mean  of  Tp 

20 

-3.78947E-06 

-1.90778E-06 

50 

1.10474E-06 

8.13508E-07 

100 

1 .68239E-06 

1,3235E-08 

200 

9  54483E-08 

3.9008E-08 

One  concern  about  the  results  in  Table  V  is  that  the  results  tended  to  consistently  get  worse 
as  we  moved  the  outer  sphere  closer  to  the  inner  sphere  of  integration.  Our  intuition  would  pre¬ 
dict  the  opposite  pattern.  The  only  apparent  explanation  seems  to  be  that  the  integral  approxi¬ 
mations  were  less  accurate  as  r  approached  R,  and  so  the  results  may  be  misleading. 

^  We  conclude  this  report  with  a  discussion  on  the  relevance  of  zero  means  in  the  functions 
Tr  and  Tr.  Often  in  physical  geodesy,  we  must  consider  the  consequences  of  the  mean  value  of 
the  geopotential  as  it  ranges  over  a  (finite)  set  of  data.  For  example,  consider  the  classical  for¬ 
mula  of  Stokes 

N  =  J  S(4i)  Ag  d0  (38) 

where  Ag  is  the  gravitational  force  field  of  the  Earth.  S  is  the  kernel,  and  N  is  the  integral  trans¬ 
formation.  Ag  has  zero  mean  over  the  global  range,  but  locally  Ag  often  fails  to  have  zero  mean. 
Furthermore,  Ag  exhibiting  zero  mean  does  not  imply  that  N  will  also  have  zero  mean.  In  prob¬ 
lems  using  (38),  it  is  sometimes  necessary  to  consider  a  global  enough  setting  so  that  the  zero 
mean  is  achieved.  Our  situation  seems  to  behave  in  a  similar  fashion.  In  the  local  setting,  the 
results  strongly  showed  that  the  estimate  T|^  and  the  actual  Tj^  do  not  have  zero  means.  In  fact, 
the  local  structure  seems  to  force  a  vast  majority  of  the  values  to  be  on  the  same  side  of  zero. 
However,  in  the  global  configuration,  we  see  a  tendency  more  toward  zero  mean  (see  Table  VI), 
but  more  sample  points  (say  1000  or  more)  would  be  needed  in  order  to  understand  how  the 
means  effect  the  results. 


VI.  SUMMARY 


1.  Geophysical  Inverse  Theory,  thus  far,  seems  a  viable  method  to  solve  the  downward 
continuation  problem,  and  map  the  geopotential  from  satellite  altitudes  to  the  Earth’s  surface. 

2.  The  method  seems  able  to  recover  the  geopotential  at  the  Earth’s  surface  with  an 
accuracy  of  a  few  percent,  with  the  absence  of  noise. 

3.  The  inverse  problem  seems  appropriately  broken  up  into  local  networks,  which  can  be 
solved  separately,  leading  to  significant  simplification  in  computing  requirements. 

4.  Further  studies  are  necessary  to  fully  understand  the  application  of  Geophysical  Inverse 
theory  to  this  problem.  Among  others  there  are: 

(a)  The  trade-off  of  size  of  the  sample  point  space  and  the  solution  space. 

(b)  The  influence  of  errors  in  the  sample  point  space,  both  position  and 
measurement  errors. 

(c)  The  benefit  of  sample  points  at  different  altitudes,  in  improving  the  solution, 
and  eliminating  possible  nonuniqueness. 

(d)  Formal  understanding  of  the  annihilator,  and  possible  tests  for  nonuniqueness. 

(e)  Applications  of  this  method  to  actual  data. 
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